Method for obtaining a phase image from an intensity image

ABSTRACT

A method of reconstructing a phase of a radiation wave field including determining a representative measure of an intensity variation of the radiation wave field on a selected surface extending globally from one end to another of the radiation wave field, determining a representative measure of the intensity of the radiation wave field on the selected surface, transforming the representative measure of the intensity variation to produce a first representation of integral transform and to apply to the first representation of integral transform a first filter and producing a first representation of modified integral transform, applying a first function of the first representation of integral transform to the first representation of modified integral transform and producing an untransformed representation, applying a correction based on the measure of intensity on the selected surface to the untransformed representation, transforming the untransformed corrected representation to produce a second representation of integral transform and applying a second filter to the second representation of integral transform, and applying a second function of the second integral transform to the second modified representation of integral transform to produce a measure of phase of the radiation wave field on the selected surface extending globally from one end to another of the radiation wave field, wherein the first and second filters are measures of a representation of a Green&#39;s function in a space of eigenfunctions of the Helmholtz equation with separable coordinates such that the first and second integral transforms are one-dimensional integral transforms in the space of the radiation wave field.

RELATED APPLICATION

This is a §371 of International Application No. PCT/FR2006/000254, withan international filing date of Feb. 3, 2006 (WO 2006/082327 A1,published Aug. 10, 2006), which is based on French Patent ApplicationNo. 05/01082, filed Feb. 3, 2005.

TECHNICAL FIELD

This disclosure relates to the field of image processing, morespecifically, to the process of obtaining an image of phase from animage of intensity.

BACKGROUND

The problem of obtaining an image of phase from an image of intensityhas been considered during the last twenty years. In <<Relationshipbetween two-dimensional intensity and phase in Fresnel Diffractionzone>>, Abramoshkin et al. (1989) provide the value of the phasegradient from intensity, the intensity gradient and Green's functions inan infinite space.

More precisely, Fornaro et al. discloses in <(<Interferometric SAR phaseunwrapping Using Green's Formulation>> using the Fast Fourier Transformsfor computing the phase from phase gradient by using Green's functions.

The formulation using the Fourier Transform is used in WO 00/26622,which describes the process of obtaining such a phase image. Thatprocess concerns the retrieval of phase of a radiation wave field bysolving the equation of the energy transfer. The rate of variation ofthe intensity of the image is determined first, orthogonally withrespect to the surface that spans the wave field (i.e., when measuringthe intensity in the two separate surfaces). This rate is then subjectto the following computation process: defining an integral transform,multiplication by a filter corresponding to the inversion of thedifferential operator, and defining the inverse integral transform. Theresult is multiplied by the function of intensity with respect to thesurface. The filters have the form based on the characteristics of theoptical system used for acquiring the intensity data, such as thenumerical aperture and spatial frequencies.

More specifically, the Fourier transform is the Fourier transform in twodimensions.

However, one of the major problems in the field of image processing isthe computation time for applications involving solution of theequations characteristic to the image.

Thus, the computation times are not sufficient to perform thecomputations for a large number of images.

Further, it should be noted that WO '622 does not mention the boundaryconditions, by supposing that the phase is zero at infinity.

SUMMARY

I provide a process of retrieval of phase of the radiation wave fieldcomprising the following steps:

-   -   compute a representative measure of the intensity variation of        the radiation wave field on a selected surface extending        globally from one end to another of the radiation wave field;    -   compute a representative measure of the intensity of the        radiation wave field on the selected surface;    -   transform the representative measure of the intensity variation        to produce a first representation of integral transform and to        apply to the first representation of integral transform a first        filter to produce a first representation of modified integral        transform;    -   apply a first function of the first representation of integral        transform to the first representation of modified integral        transform to produce an untransformed representation;    -   apply a correction based on the measure of intensity on the        selected surface to the untransformed representation;    -   transform the untransformed corrected representation to produce        a second representation of integral transform and to apply to        the second representation of integral transform a second filter;    -   apply a second function of the second integral transform to the        second modified representation of integral transform to produce        a measure of phase of the radiation wave field on the selected        surface extending globally from one end to another of the        radiation wave field,    -   wherein the first and second filters are measures of a        representation of a Green's function in a space of        eigenfunctions of the Helmholtz equation with separable        coordinates such that the first and second integral transforms        are one-dimensional integral transforms in the space of the        radiation wave field.

The first function may use a trigonometric function of the cosine type,and the second function uses a trigonometric function of the sine type.

By the same token, it is possible that the first function uses atrigonometric function of the sine type, and the second function uses atrigonometric function of the cosine type.

As an advantage, the first and second integral transforms are producedby using a Fourier transform, for example, a fast transform. Accordingto one alternative, the first and second filters are sensibly the same.The surfaces may be sensibly planar.

I also provide a computer program, possibly saved on a support, toexecute the different steps of the process.

I further provide an apparatus of reconstruction of phase of a radiationwave field comprising:

-   -   a means for acquisition of a representative measure of the        intensity variation of the radiation wave field on a selected        surface extending globally from one end to another of the        radiation wave field;    -   a means for acquisition of a representative measure of the        intensity of the radiation wave field on the said selected        surface;    -   a processing means to sequentially perform the following steps:        -   (i) transform the representative measure of the intensity            variation to produce a first representation of integral            transform and to apply to the first representation of            integral transform a first filter to produce a first            representation of modified integral transform;        -   (ii) apply a first function of the first representation of            integral transform to the first representation of modified            integral transform to produce an untransformed            representation;        -   (iii) apply a correction based on the measure of intensity            on the selected surface to the untransformed representation;        -   (iv) transform the untransformed corrected representation to            produce a second representation of integral transform and to            apply to the second representation of integral transform a            second filter;        -   (v) apply a second function of the second integral transform            to the second modified representation of integral transform            to produce a measure of phase of the radiation wave field on            the selected surface extending globally from one end to            another of the radiation wave field,        -   wherein the first and second filters are measures of a            representation of a Green's function in a space of            eigenfunctions of the Helmholtz equation with separable            coordinates such that the first and second integral            transforms are one-dimensional integral transforms in the            space of the radiation wave field.

The means of acquisition of a representative measure of the intensityvariation of the radiation wave field may comprise at least one mobileplatform. The means of acquisition of a representative measure of theintensity variation of the radiation wave field may also comprise atleast a lens with a variable length of focus.

BRIEF DESCRIPTION OF THE DRAWINGS

My disclosure is better understood due to the description made hereafterfor purely explanatory purposes of a representative example, withreference to the annexed figures where:

FIG. 1 represents image data;

FIG. 2 shows a comparison between the number of complex multiplicationsin my methods and in WO 00/26622 for a processed image with dimensionsfrom 64×64 pixels to 128×128 pixels;

FIG. 3 shows a comparison between the number of complex multiplicationsin my methods and in WO 00/26622 for a processed image with dimensionsfrom 1024×1024 pixels to 2047×2047 pixels;

FIG. 4 shows the data from initial images at the start of the algorithm;

FIGS. 5A, 5B and 5C, which should be read in a successive manner,represent one exemplary implementation mode; and

FIGS. 6-9 represent devices for implementation of my methods.

My disclosure is described in detail below in the Cartesian coordinatesystem and with a rectangular image. It is understood that one skilledin the art is capable, in reading this description, to adapt theinformation contained herein for the other types of images and/or theother types of coordinates.

DETAILED DESCRIPTION

I provide a method of reconstruction of the function φ({right arrow over(r)}) of phase in a point {right arrow over (r)}=(x, y) of a 2D planefrom samples of the first-order differences in points (x, y)=(nΔx, mΔy).We compute the phase only in the points for which x=nΔx and y=mΔy, wherethe lines of the two-dimensional (2D) grid intersect, i.e., in thecenters of points where elements (pixels) of a two-dimensional sensor ofthe luminous energy are situated, for example, of a CCD camera. Thephase differences of the first order in the points (x, y)=(nΔx, mΔy) areobtained by computation of the gradient of the phase from the derivativeof the intensity ∂I/∂z along the optical axis z. This computation isdetailed later. To facilitate the transition from the discrete notationto the continuous notation, the description of terms used in thedescription of the process is presented in the Table I.

Fornaro et al. used the first Green's identity for evaluation of thephase from the differences of the first order presented in thecontinuous form as a gradient vector {right arrow over (F)}({right arrowover (r)}). The phase is computed as a sum of the integral over thesupport region S of the image of phase and over the boundary C of thisregion:φ({right arrow over (r)}′)=∫_(S) {right arrow over (F)}({right arrowover (r)})·∇g({right arrow over (r)},{right arrow over (r)}′)d{rightarrow over (r)}−

φ( {right arrow over (r)})[{circumflex over (n)} _(S) ·∇g({right arrowover (r)},{right arrow over (r)}′)]d{right arrow over (c)}  Equation 1where {circumflex over (n)}_(S) is a unity vector perpendicular to C andoriented exterior to S as shown in FIG. 1, and where g({right arrow over(r)}, {right arrow over (r)}′) is a Green's function. In Equation 1, thephase φ enters both in the left-hand side and in the integral on thecontour C. This double introduction of the same value in the two partsof Equation 1 does not obtain the phase if it is not zero or known apriori on the contour C.

Therefore, the integration of the phase on the contour C means that thephase should either be known, be equal to zero, or be pre-calculated byanother optical device, for example, the Shack-Hartmann arrangement oflenses, prior to reconstructing the phase inside the surface S. Theintegration on the contour in Equation 1 can also be eliminated if thechosen Green's function satisfies the Neumann boundary conditions on thecontour C of the surface S:{circumflex over (n)} _(S) ·∇g({right arrow over (r)},{right arrow over(r)}′)=0  Equation 2

Equation 2 is valid for the Green's function, and therefore the phasereconstructed according to the Equation 1 satisfies Neumann boundaryconditions. The integration on the contour C can also be eliminated ifthe chosen Green's function satisfies, on the contour C, Dirichletboundary conditions and if it is known that the phase is zero on C, thatis written as:g({right arrow over (r)},{right arrow over (r)}′)=0,{right arrow over(r)}εC,{right arrow over (r)}′εSφ({right arrow over (r)})=0,{right arrow over (r)}εC.  Equation 3

As previously, Equation 3 is valid for the Green's function and for thephase on the contour C, therefore the phase computed according to (1)using the Green's function from Equation 3 satisfies Dirichlet boundaryconditions.

Further, it is known that the Green's function satisfies the equation:∇² g({right arrow over (r)},{right arrow over (r)}′)=−δ({right arrowover (r)},{right arrow over (r)}′)  Equation 4where δ({right arrow over (r)},{right arrow over (r)}′) is the Diracfunction.

To satisfy preceding equations, a Green's function may be chosen to beexpressed in the form of a series in eigenfunctions of the Helmholtzequation as follows. One skilled in the art will note that itautomatically satisfies the boundary conditions (2) or (3) according tothe chosen form, and also the identity (4).

For a rectangular support region S, these special functions are:U _(m,n)({right arrow over (r)})=A _(m,n) cos(ξ_(n) x)cos(η_(m) y),m,n=0, . . . , ∞  Equation 5for the Neumann boundary conditions, andU _(m,n)({right arrow over (r)})=A _(m,n) sin(ξ_(n) x)sin(η_(m) y),m,n=0, . . . , ∞  Equation 6for the Dirichlet boundary conditions.

For a support S of dimensions a×b, the constant A_(m,n) in 5 and 6 isequal to √{square root over (2)}a^(−1/2)b^(−1/2) if m or n are zero and2a^(−1/2)b^(−1/2) otherwise, where a is the length, b is the height ofthe rectangular region of support S as in FIG. 1.

The eigenvalues associated with the functions (5) and (6) are also used,which are:

$\begin{matrix}{{\sigma_{m,n}^{2} = {\xi_{n}^{2} + \eta_{m}^{2}}},{\xi_{n} = \frac{n\;\pi}{a}},{\eta_{m} = \frac{m\;\pi}{b}},m,{n = 0},\ldots\mspace{11mu},{\infty.}} & {{Equation}\mspace{20mu} 7}\end{matrix}$

The Green's function for the domain S is written as an infinite seriesin eigenfunctions of the Helmholtz equation:

$\begin{matrix}{{g\left( {\overset{\rightarrow}{r},{\overset{\rightarrow}{r}}^{\prime}} \right)} = {\sum\limits_{m,{n = 0},{\sigma_{m,n}^{2} \neq 0}}^{\infty}{\sigma_{m,n}^{- 2}{U_{m,n}\left( \overset{\rightarrow}{r} \right)}{{U_{m,n}\left( {\overset{\rightarrow}{r}}^{\prime} \right)}.}}}} & {{Equation}\mspace{20mu} 8}\end{matrix}$

The derivation of the general expression (8) and the verification of theboundary conditions (2) and (3) are readily realizable by one skilled inthe art for such a function.

Thus, if the Green's function (8) is used as g({right arrow over(r)},{right arrow over (r)}′) in the Equation 1, the integration on thecontour in the Equation 1 can be suppressed and the phase can beexpressed in the compact form:φ({right arrow over (r)}′)=∫_(S) {right arrow over (F)}({right arrowover (r)})·∇g({right arrow over (r)},{right arrow over(r)}′)dS  Equation 9in which it is noted that the gradient is computed with respect to thevariable {right arrow over (r)}, that is ∇=∇_({right arrow over (r)}).

One aspect is to compute numerically the Equation 9 using the definedGreen's functions to obtain the required phase, in preserving theacceptable computation time. The solution is obtained by the possibilityto use a one-dimensional Fast Fourier Transform (1DFFT).

The use of this transform in one dimension is not straightforward forone skilled in the art taking into consideration general equationsabove. I provide below the various steps for comprehension of the factthat it is possible to use such a transform in one dimension.

Several preliminary transformations are necessary. To facilitatecalculations, I point out various intermediate stages of calculationwhich should not be regarded as essential components, but simply themeans of obtaining a form of equation with which it is possible toimplement the transforms in an adequate dimension.

For more clarity, the Neumann boundary conditions are considered first,with γ=a/b and the Green's function in the expression of the Equation 8with special functions (5):

$\begin{matrix}{{g\left( {\overset{\rightarrow}{r},{\overset{\rightarrow}{r}}^{\prime}} \right)} = {\sum\limits_{{m = 0},{{m^{2} + n^{2}} \neq 0}}^{\infty}{A_{m,n}^{2}{\cos\left( {\eta_{m}y} \right)}{\cos\left( {\eta_{m}y^{\prime}} \right)} \times \mspace{326mu}{\sum\limits_{n = 0}^{\infty}{\frac{{\cos\left( {\xi_{n}x} \right)}{\cos\left( {\xi_{n}x^{\prime}} \right)}}{{\left( {\pi/a} \right)n^{2}} + {\left( {\pi/b} \right)m^{2}}}.}}}}} & {{Equation}\mspace{20mu} 10}\end{matrix}$

The cases n=0, m≠0 and n≠0, m=0 are separated to obtain:

$\begin{matrix}{{g\left( {\overset{\rightarrow}{r},{\overset{\rightarrow}{r}}^{\prime}} \right)} = {{2\pi^{- 2}\gamma{\sum\limits_{n = 1}^{\infty}{n^{- 2}{\cos\left( {\xi_{n}x} \right)}{\cos\left( {\xi_{n}x^{\prime}} \right)}}}} + {2\pi^{- 2}\gamma^{- 1}{\sum\limits_{m = 1}^{\infty}{m^{- 2}{\cos\left( {\eta_{m}y} \right)}{\cos\left( {\eta_{m}y^{\prime}} \right)}}}} + {4a^{- 1}b^{- 1}{\sum\limits_{m = 1}^{\infty}{{\cos\left( {\eta_{m}y} \right)}{\cos\left( {\eta_{m}y^{\prime}} \right)}{\sum\limits_{n = 1}^{\infty}{\frac{{\cos\left( {\xi_{n}x} \right)}{\cos\left( {\xi_{n}x^{\prime}} \right)}}{\left( {n/a} \right)^{2} + \left( {m/b} \right)^{2}}.}}}}}}} & {{Equation}\mspace{20mu} 11}\end{matrix}$

However, it is known that:

$\begin{matrix}{{{\sum\limits_{n = 1}^{\infty}\frac{\cos({nu})}{n^{2} + \alpha^{2}}} = {{\frac{\pi}{2\alpha}\frac{\cosh\left\lbrack {\alpha\left( {\pi - u} \right)} \right\rbrack}{\sinh({\alpha\pi})}} - \frac{1}{2\alpha^{2}}}},{0 \leq u \leq {2\pi}},} & {{Equation}\mspace{20mu} 12}\end{matrix}$

Therefore,

$\begin{matrix}{{{a^{- 1}b^{- 1}{\sum\limits_{n = 1}^{\infty}\frac{{\cos\left( {\xi_{n}x} \right)}{\cos\left( {\xi_{n}x^{\prime}} \right)}}{\left( {n/a} \right)^{2} + \left( {m/b} \right)^{2}}}} = {{\frac{\pi}{2m}{f_{m}^{\gamma}\left( {x,x^{\prime}} \right)}} - \frac{1}{2\gamma\; m^{2}}}},\mspace{20mu}{m \neq 0},\mspace{20mu}{with}} & {{Equation}\mspace{20mu} 13} \\{{f_{m}^{\gamma}\left( {x,x^{\prime}} \right)} = {\frac{1}{\sinh\left( {\gamma\; m\;\pi} \right)}\left\{ \begin{matrix}{{{\cosh\left\lbrack {\eta_{m}\left( {a - x} \right)} \right\rbrack}{\cosh\left( {\eta_{m}x^{\prime}} \right)}},} & {{{if}\mspace{14mu} x} > x^{\prime}} \\{{{\cosh\left\lbrack {\eta_{m}\left( {a - x^{\prime}} \right)} \right\rbrack}{\cosh\left( {\eta_{m}x^{\prime}} \right)}},} & {{{if}\mspace{14mu} x} < {x^{\prime}.}}\end{matrix} \right.}} & {{Equation}\mspace{20mu} 14}\end{matrix}$

However,

$\begin{matrix}{{\sum\limits_{n = 1}^{\infty}{n^{- 2}{\cos({nu})}{\cos({nv})}}} = {\frac{1}{12}\left\{ \begin{matrix}{{{3u^{2}} + {3\left( {v - \pi} \right)^{2}} - \pi^{2}},} & {{{if}\mspace{14mu} 0} \leq u \leq v} \\{{{3v^{2}} + {3\left( {u - \pi} \right)^{2}} - \pi^{2}},} & {{{if}\mspace{14mu} 0} \leq u \leq {v.}}\end{matrix} \right.}} & {{Equation}\mspace{20mu} 15}\end{matrix}$

Therefore, the Green's function is given according to this aspect by:

$\begin{matrix}{{g\left( {\overset{\rightarrow}{r},{\overset{\rightarrow}{r}}^{\prime}} \right)} = {{2\pi^{- 1}{\sum\limits_{m = 1}^{\infty}{m^{- 1}{\cos\left( {\eta_{m}y} \right)}{\cos\left( {\eta_{m}y^{\prime}} \right)}{f_{m}^{\gamma}\left( {x,x^{\prime}} \right)}}}} + {\frac{\gamma}{6a^{2}}\left\{ \begin{matrix}{{{3x^{2}} + {3\left( {x^{\prime} - a} \right)^{2}} - a^{2}},} & {{{if}\mspace{14mu} x} \leq x^{\prime}} \\{{{3\left( x^{\prime} \right)^{2}} + {3\left( {x - a} \right)^{2}} - a^{2}},} & {{{if}\mspace{14mu} x^{\prime}} < {x.}}\end{matrix} \right.}}} & {{Equation}\mspace{20mu} 16}\end{matrix}$

However, it is also known that:

$\begin{matrix}{{g\left( {\overset{\rightarrow}{r},{\overset{\rightarrow}{r}}^{\prime}} \right)} = {{2\pi^{- 1}{\sum\limits_{m = 1}^{\infty}{n^{- 1}{\cos\left( {\xi_{n}x} \right)}{\cos\left( {\xi_{n}x^{\prime}} \right)}{f_{m}^{\overset{\_}{\gamma}}\left( {y,y^{\prime}} \right)}}}} + {\frac{\overset{\_}{\gamma}}{6b^{2}}\left\{ \begin{matrix}{{{3y^{2}} + {3\left( {y^{\prime} - b} \right)^{2}} - b^{2}},} & {{{if}\mspace{14mu} y} \leq y^{\prime}} \\{{{3\left( y^{\prime} \right)^{2}} + {3\left( {y - b} \right)^{2}} - b^{2}},} & {{{{if}\mspace{14mu} y^{\prime}} < y},}\end{matrix} \right.}}} & {{Equation}\mspace{20mu} 17}\end{matrix}$where γ=γ⁻¹ with:

$\begin{matrix}{{f_{n}^{\overset{\_}{\gamma}}\left( {y,y^{\prime}} \right)} = {\frac{1}{\sinh\left( {\overset{\_}{\gamma}\; n\;\pi} \right)}\left\{ \begin{matrix}{{{\cosh\left\lbrack {\xi_{n}\left( {b - y} \right)} \right\rbrack}{\cosh\left( {\xi_{n}y^{\prime}} \right)}},} & {{{if}\mspace{14mu} y} > y^{\prime}} \\{{{\cosh\left\lbrack {\xi_{n}\left( {b - y^{\prime}} \right)} \right\rbrack}{\cosh\left( {\xi_{n}y^{\prime}} \right)}},} & {{{if}\mspace{14mu} y} < {y^{\prime}.}}\end{matrix} \right.}} & {{Equation}\mspace{20mu} 18}\end{matrix}$

By taking the derivatives of these expressions of the Green's functionwith respect to x and y, it follows that:

$\begin{matrix}{{{g_{y}\left( {\overset{\rightarrow}{r},{\overset{\rightarrow}{r}}^{\prime}} \right)} = {{- 2}b^{- 1}{\sum\limits_{m = 1}^{\infty}{{\sin\left( {\eta_{m}y} \right)}{\cos\left( {\eta_{m}y^{\prime}} \right)}{f_{m}^{\gamma}\left( {x,x^{\prime}} \right)}}}}}{and}} & {{Equation}\mspace{20mu} 19} \\{{g_{x}\left( {\overset{\rightarrow}{r},{\overset{\rightarrow}{r}}^{\prime}} \right)} = {{- 2}a^{- 1}{\sum\limits_{n = 1}^{\infty}{{\sin\left( {\xi_{n}x} \right)}{\cos\left( {\xi_{n}x^{\prime}} \right)}{f_{n}^{\overset{\_}{\gamma}}\left( {y,y^{\prime}} \right)}}}}} & {{Equation}\mspace{20mu} 20}\end{matrix}$and in inserting equation (9) it follows that:

$\begin{matrix}{{\phi\left( {x^{\prime},y^{\prime}} \right)} = {{2b^{- 1}{\sum\limits_{m = 1}^{M^{\prime} - 1}{{\cos\left( {\eta_{m}y^{\prime}} \right)}{\int_{0}^{a}{{f_{m}^{\gamma}\left( {x,x^{\prime}} \right)}{D_{x}\left( \eta_{m} \right)}{\mathbb{d}x}}}}}} + {2a^{- 1}{\sum\limits_{n = 1}^{N^{\prime} - 1}{{\cos\left( {\xi_{x}x^{\prime}} \right)}{\int_{0}^{b}{{f_{n}^{\overset{\_}{\gamma}}\left( {y,y^{\prime}} \right)}{D_{y}\left( \xi_{n} \right)}{\mathbb{d}y}}}}}}}} & {{Equation}\mspace{20mu} 21}\end{matrix}$where M′ and N′ are integers. They can be chosen to attain the requirednumerical precision. In particular, I choose these two numbers as powersof two, greater than M+1 and N+1, respectively, where M and N representthe dimensions of the image to process.

Due to the presence of the factors sin(η_(m)y) and sin(ξ_(n)x) inequations (19) and (20), the integration [with respect to y in the firstsum of the right-hand side of Equation 21 to obtain D_(x)(η_(m)) andwith respect to x in the second sum of the right-hand side to obtainD_(y)(ξ_(n))] allows one to compute the sequences D_(x)(η_(m)) andD_(y)(ξ_(n)) as the imaginary parts (taken with the negative sign) ofthe Fourier Transform (FT) of the derivatives (with respect to y andwith respect to x, respectively) of the phase.

These two sequences can be computed numerically by using the algorithmof the Discrete Fourier Transform, by algorithm of the Fast FourierTransform (1DFFT or FFT 1D), as applied to the phase differences of thefirst order F_(y) and F_(x) respectively:D _(x)(η_(m))=−

{FFT└F _(y)(n,m),n=const┘}_(η) _(m)   Equation 22D _(y)(ξ_(n))=−

{FFT[F _(x)(n,m),m=const]}_(ξ) _(n)   Equation 23

It is noted that n=const means that the first index (i.e., n) is fixedin F_(y)(n, m) whereas the FFT 1D of the resultant sequence, in whichthe second index (i.e., m) varies, is computed. Inversely, m=const meansthat the second index (i.e., m) is fixed in F_(x)(n, m) whereas the FFT1D of the resultant sequence, in which the first index (i.e., n) varies,is computed. The indices η_(m) and ξ_(n) indicate the discretefrequencies for which the 1D FFT should be evaluated.

In addition, taking into consideration calculation above, in case ofDirichlet boundary conditions, the phase is given in the same way by:

$\begin{matrix}{{\phi\left( {x^{\prime},y^{\prime}} \right)} = {{2b^{- 1}{\sum\limits_{m = 1}^{M^{\prime} - 1}{{\sin\left( {\eta_{m}y^{\prime}} \right)}{\int_{0}^{a}{{f_{m}^{\gamma}\left( {x,x^{\prime}} \right)}{D_{x}\left( \eta_{m} \right)}{\mathbb{d}x}}}}}} + {2a^{- 1}{\sum\limits_{n = 1}^{N^{\prime} - 1}{{\sin\left( {\xi_{n}x^{\prime}} \right)}{\int_{0}^{b}{{f_{m}^{\overset{\_}{\gamma}}\left( {y,y^{\prime}} \right)}{D_{y}\left( \xi_{n} \right)}{\mathbb{d}y}}}}}}}} & {{Equation}\mspace{20mu} 24}\end{matrix}$with the real parts of the Fourier transforms in one dimension:D _(x)(η_(m))=

{FFT└F _(y)(n,m), n=const┘}_(η) _(m)D _(y)(ξ_(n))=

{FFT[F _(x)(n,m), m=const]}_(ξ) _(n) .  Equation 25Computation From The Intensity

Again with reference to Abramoshkin et al., the gradient of the phasefrom the Green's function is known:

$\begin{matrix}{{{\nabla{\phi\left( {\overset{\rightarrow}{r}}^{\prime} \right)}} = {\int_{S}{\frac{\partial{I\left( \overset{\rightarrow}{r} \right)}}{\partial z}{\nabla{g\left( {\overset{\rightarrow}{r},{\overset{\rightarrow}{r}}^{\prime}} \right)}}{\mathbb{d}\overset{\rightarrow}{r}}}}},} & {{Equation}\mspace{20mu} 26}\end{matrix}$where the Green's function used is the Green's function for the infinitespace:

$\begin{matrix}{{g\left( {\overset{\rightarrow}{r},{\overset{\rightarrow}{r}}^{\prime}} \right)} = {{- \frac{1}{2\pi}}\ln{{{\overset{\rightarrow}{r} - {\overset{\rightarrow}{r}}^{\prime}}}.}}} & {{Equation}\mspace{20mu} 27}\end{matrix}$

Therefore, by re-using Equation 1,

$\begin{matrix}{{\phi\left( {\overset{\rightarrow}{r}}^{\prime} \right)} = {{{- k}{\int_{S}{\left\lbrack {I\left( {\overset{\rightarrow}{r};z} \right)} \right\rbrack^{- 1}{\left\{ {\int_{S}{\frac{\partial{I\left( {\overset{\rightarrow}{r}}^{''} \right)}}{\partial z}{\nabla{g\left( {{\overset{\rightarrow}{r}}^{''},\overset{\rightarrow}{r}} \right)}}{\mathbb{d}{\overset{\rightarrow}{r}}^{''}}}} \right\} \cdot {\nabla{g\left( {\overset{\rightarrow}{r},{\overset{\rightarrow}{r}}^{\prime}} \right)}}}{\mathbb{d}\overset{\rightarrow}{r}}}}} + {\oint_{C}{{{\phi\left( \overset{\rightarrow}{r} \right)}\left\lbrack {{\hat{n}}_{S} \cdot {\nabla{g\left( {\overset{\rightarrow}{r},{\overset{\rightarrow}{r}}^{\prime}} \right)}}} \right\rbrack}{{\mathbb{d}\overset{\rightarrow}{c}}.}}}}} & {{Equation}\mspace{20mu} 28}\end{matrix}$

By using the Green's function (that cancel the second integral on thecontour), it follows that:

$\begin{matrix}{{\phi\left( {\overset{\rightarrow}{r}}^{\prime} \right)} = {{- k}{\int_{S}{\left\lbrack {I\left( {\overset{\rightarrow}{r};z} \right)} \right\rbrack^{- 1}{\left\{ {\int_{S}{\frac{\partial{I\left( {\overset{\rightarrow}{r}}^{''} \right)}}{\partial z}{\nabla{g\left( {{\overset{\rightarrow}{r}}^{''},\overset{\rightarrow}{r}} \right)}}{\mathbb{d}{\overset{\rightarrow}{r}}^{''}}}} \right\} \cdot {\nabla{g\left( {\overset{\rightarrow}{r},{\overset{\rightarrow}{r}}^{\prime}} \right)}}}{\mathbb{d}\overset{\rightarrow}{r}}}}}} & {{Equation}\mspace{20mu} 29}\end{matrix}$with, for Neumann conditions:

$\begin{matrix}{{{\phi_{y^{\prime}}\left( {x^{\prime},y^{\prime}} \right)} = {2b^{- 1}{\sum\limits_{m = 1}^{M^{\prime} - 1}{{\cos\left( {\eta_{m}y^{\prime}} \right)}{\int_{0}^{a}{{f_{m}^{\gamma}\left( {x,x^{\prime}} \right)}{P_{x}\left( \eta_{m} \right)}{\mathbb{d}x}}}}}}},{{\phi_{x^{\prime}}\left( {x^{\prime},y^{\prime}} \right)} = {2a^{- 1}{\sum\limits_{n = 1}^{N^{\prime} - 1}{{\cos\left( {\xi_{n}x^{\prime}} \right)}{\int_{0}^{b}{{f_{m}^{\overset{\_}{\gamma}}\left( {y,y^{\prime}} \right)}{P_{y}\left( \xi_{n} \right)}{\mathbb{d}y}}}}}}},\mspace{20mu}{with}} & {{Equation}\mspace{20mu} 30} \\{\mspace{79mu}{{{P_{x}\left( \eta_{m} \right)} = {{- {??}}\left\{ {{FFT}\left\lbrack {{\frac{\partial I}{\partial z}\left( {n,m} \right)},{n = {const}}} \right\rbrack} \right\}_{\eta_{m}}}}\mspace{20mu}{{P_{y}\left( \xi_{n} \right)} = {{- {??}}{\left\{ {{FFT}\left\lbrack {{\frac{\partial I}{\partial z}\left( {n,m} \right)},{m = {const}}} \right\rbrack} \right\}_{\xi_{n}}.}}}}} & {{Equation}\mspace{20mu} 31}\end{matrix}$

One also obtains the equivalent functions for the Dirichlet boundaryconditions.

Equation 28 is substantially that used in WO '622. In WO '662, thesecond term is cancelled by supposing that the phase is zero on thecontour, whereas my phase can be nonzero. The calculation isnevertheless realizable by Neummann and Dirichlet boundary conditionscancelling this term. In addition, in WO '622, Equation 28 is expressedin the Fourier domain and calculated by Fourier transforms of Fourier intwo dimensions. On the contrary, I carry out this calculation of phaseby using the Fourier transforms in one dimension thanks to the use ofspecific Green's functions which substantially decreases the number ofcomplex operations, and thus the computing time of the algorithmsassociated with my method.

As illustrated in FIGS. 2 and 3, the number of complex multiplicationsis calculated between my solution by the 1D Fourier transform, and thesolution of WO '622 depending on the size of the processed image. It isclear that my solution saves considerably the image processing time forthe phase reconstruction.

Processing Algorithm for Phase Reconstruction

I now provide an example. For this, the following notations are used:

Abbreviations:

-   -   BC—Boundary Conditions;    -   1D FFT—One-dimensional direct Fourier Transform;    -   1D IFFT—One-dimensional inverse Fourier Transform;    -   [ . . . ]—Real part of the complex number in brackets;    -   [ . . . ]—Imaginary part of the complex number in brackets.        The following parameters of the initial images:    -   M×N, where M is the number of rows in the image, indexed by i=1,        . . . , M; and N is the number of columns, indexed by j=1, . . .        , N.        Constants:    -   a=N′=2^(i)≧N+1;    -   b=M′=2^(j)≧M+1;    -   H_(max)=225 (for the 32-bit computer);    -   γ=a/b, Δν_(y)=π/b, M_(max)=min{M′, H_(max)/γ};    -   β=b/a, Δν_(x)=π/a, N_(max)=min{N′, H_(max)/β}.

FIG. 4 shows the first step of the algorithm in which, from an image,one computes ∂I/∂z by difference between intensity images in the planesz₁ and z₂. This difference is of course computed pixel by pixel in theimages I₂(m, n) and I₁(m, n). It will be understood that, in case ofarbitrary coordinates, this measure of ∂I/∂z corresponds in fact to anymeasure of the intensity variation over a surface.

In a general way, the algorithm gives the possibility, each time whenthis is necessary, of choosing between Neumann and Dirichlet boundaryconditions.

FIGS. 5 a, 5 b and 5 c show an implementation of the equations describedabove. It is also noted that the discrete nature of the real data in theimages (values in x and y being obtained from pixel data), allows adiscrete calculation of the integrals given in the equations and anformulation of calculations in the form of finite discrete sums. Themodes of calculations of these discrete sums are given in the algorithmas an example and can of course vary according to choices withoutlimiting the scope of the appended claims.

Illustrated in FIG. 5 a, the oblique hatchings of the initial imagedenote the FFT window used for calculations. This window exceeds theuseful image and can be filled by 0, according to the known method ofthe “zero padding”, or “zero-filling”, or by any other well known methodin the field of the image processing.

The vertical lines of FIG. 1 indicate the direction of the 1D FFT inthis illustrated part of the algorithm. In this mode of implementation,the FFT is calculated by fixing a column while varying the row index.

In addition, it can be noted from FIG. 5 a the introduction of theparameter M_(max) used preferably according to this implementation modeto avoid divergence of the terms in hyperbolic sine and a specificcalculation of the coefficient c_(y)(m,n) if the index representing thecoordinates exceeds this threshold value.

One will recognize, by Operation A, obtaining the phase derivative withrespect to y as described in the first part of Equation 30. It is notedthat Operation A requires the operation of division by a measure ofintensity which can make the solution diverge if this intensity is zero.In this case, it is possible to carry out a determination of the zerosof the intensity beforehand to correct the divergence and/or to assign aminimum value to the intensity in the event of zero intensity in aconsidered point of space.

According to Equation 22, the method of calculation applies a 1D Fouriertransform to the phase gradient with respect to y to obtain the firstpart of the phase as in Equation 21. This first term φ₂(m, n) isobtained as a result of the Operation B of FIG. 5 c.

In the same manner (not shown), the phase derivative with respect to xis obtained as described in the first part of Equation 30, and accordingto Equation 23, one applies then the method of calculation using the 1DFourier transform applied to the phase gradient with respect to x toobtain the second part of the phase φ₁(m, n), as in Equation 21.

The reconstructed phase is computed by φ(m, n)=φ₁(m, n)+φ₂(m, n).

In a general way, illustrated in FIGS. 5A, 5B and 5C, the algorithm thuscomprises, after the calculation of a measure of intensity variation ona surface (∂I/∂z in the illustrated realization mode), the applicationof a first integral transformation, preferably in the form of a FastFourier Transform, to obtain the variable named r_(x)(m, n) depending onthe used boundary conditions. This variable is then multiplied (in theFourier space or, more generally, in the projection space) by a filtercorresponding to the functions ƒ_(n) ^(γ) and ƒ_(m) ^(γ), possiblymodified by the calculation constraints of the type M_(max).

In accordance with Equations 19 and 20, these filters ƒ_(n) ^(γ) andƒ_(m) ^(γ) are representations of a Green's function in a space ofeigenfunctions of the Helmholtz equation, for example the functionsassociated with Equations 5 and 6.

Inverse integral transforms are then applied, for example, InverseFourier Transforms, by taking the real or imaginary parts according toboundary conditions.

For example, as Equation 26 indicates, the division by a measure ofintensity, for example, that taken on a computation plane of theintensity variation, makes it possible to obtain the gradient of thephase with respect to one of the coordinates.

In conformity with Equations 21, 22, and 23, a one-dimensional integraltransform, for example, of Fourier type, is then applied, a secondfilter to the obtained variable [D_(x)(m, n) in the proposed realizationmode] is applied, and an inverse integral transform is applied, toobtain substantially one of the terms of the expression defining therequired phase.

It is understood that using the example of implementation of thealgorithm given in the figures and described above, it is to makemodifications of implementation, for example, according to computationpower of the used computers or according to specific characteristics ofthe image.

Regarding the description made above of the equations necessary for theimplementation of the algorithm and of the processing algorithm itself,it will be understood that the described realization mode is only anon-restrictive example.

In particular, in the case of a nonrectangular support, it is possibleto determine the eigenfunctions of the Helmholtz equation for a specialform of support. As an indication only, the form of these functions incase of a circular support is given by:U _(m,n)(ρ,φ)=[cos(mφ)+sin(mφ)]J _(m)(πα_(mn) ρ/a), m, n=0, . . . ,∞  Equation 5where J_(m) is the Bessel function or order m with J′_(m)(πα_(mn))=0, m,n=0, . . . , ∞, i.e. πα_(mn) is the n-th zero of the derivative ofBessel function of order m. Equation 8 allows then to obtain theassociated Green's function.

The use of the Green's function associated with this eigenfunction isparticularly useful if the phase reconstruction in a circular area givesthe sufficient information about the topography of the object ofcircular geometry, such as for example in case of visualization andmeasurement of the front face of an optical fiber.

One skilled in the art will be able to determine these eigenfunctions inthe case of triangular or semi-infinite supports.

In addition, I describe an implementation mode utilizing one-dimensionaldiscrete Fourier transforms, but one skilled in the art will be able toimplement one-dimensional discrete Bessel transforms, as a curvilinearanalogue of the Fourier transforms. It is also possible to use any typeof integral transformation as a projection of a function on a functionspace, for example, determined by the type of function used in equations5 or 5′.

In any case, if the eigenfunctions of Helmholtz equation are expressedin separable coordinates, in particular if the support of the image isnot distorted too much, my methods can be used and integraltransformations can be carried out in only one dimension.

For very distorted supports, conformal mapping transforms can be used toreduce these supports to simpler supports, or to define this support asa composition of simpler supports.

Moreover, it is understood that the implementation mode described for animage is also applicable for any type of radiation wave field in thewhole luminous spectrum. To apply my methods at the various spectra, itis possible to modify the acquisition setups as described in thefollowing section devoted to the setups associated with the invention.

Measurement Setups

In accordance with the example in FIG. 4, it is advisable to use devicesof measurement of the intensity in distinct planes. These measurementsare carried out in planes perpendicular to the optical axis.

To do this, as illustrated in FIG. 6, the acquisition of an image of atransparent sample is carried out. Measurements can be taken on the twoperpendicular planes, the position of the sample can be modifiedcompared to a fixed optical device or to move the optical device.

In FIG. 6, for a transparent sample, the device comprises an imagesensor 61, for example, in form of a CCD sensor, an optical projectionsystem 62 to project the light coming from sample 63 towards the sensor.The light is emitted by an optical illuminating system 64. The sample ismobile along the Z axis due to a displacement device (not shown).

As shown in FIG. 7, the sample is fixed at a motionless platform 75. Theoptical device still includes a sensor 71, a projection system 72 and anillumination system 73, and it includes also a movable platform 76 towhich elements 71, 72 and 73 are attached. This platform is set inmotion by a displacement system 77 in form of, for example, a motor or apiezoelectric element.

As shown in FIG. 8, I include a measurement device using one or morelenses with variable focal distances. In this case, a lens with avariable focal distance 88 is placed for example at the position of theprojection system and, according to an alternative, a lens with variablefocal distance 89 is placed at the position of the illumination systemfor measurements of a fixed transparent sample. In this case, thedisplacement in intensity measurements corresponds to the difference infocal distances, Δz=F₂−F₁. In case of two lenses, the modifications ofthe focal distances will be done simultaneously for the two lenses. Thissystem is also realizable with only one lens with a variable focaldistance.

As shown in FIG. 9, for the diffusing and non-transparent sample, onecan use a beam splitter 90 intended to direct a luminous beam comingfrom the illuminating system 94. One can then use the set of theelements previously described such as the sensor, the projection system,the whole unit being either fixed or mobile, as well as one or morelenses with variable focal distances, as previously described.

The image sensor is connected to an image processing device including asoftware for the implementation of the algorithm.

The images of phase obtained by this software can then be visualized ona computer screen.

TABLE I Term in the continuous Description of the term notation in thecontinuous notation Relation to the discrete notation Relation to thegeometry of a CCD type sensor r = (x, y) and Vectors from the origin tothe points Grid of points with coordinates (m, n) for each Coordinatesof the centre of an element on the r′ = (x′, y′) in the two-dimensionalspace integer number m and n. CCD sensor S Support region for the phaseGrid of points with coordinates (m, n) where All the active elements ofthe sensor m = 1, . . . , M and n = 1, . . . , N. x = a Boundary line tothe phase support The column of the grid with the column index Thecolumn of the passive elements of the region S, in the y direction n =N + 1 sensor situated to the right of all the active elements of thesensor. y = b Boundary line to the phase support The row of the gridwith the row index The row of the passive elements of the sensor regionS, in the x direction m = M + 1 situated beneath all the active elementsof the sensor C Boundary to the phase support Two rows of the grid ofpoints with coordinates All the passive elements of the sensor (or itsregion S (0, n) and (M + 1, n) where n = 1, . . . , N, boundary) thatenclose the rectangular area of as well as two columns with coordinates(m, 0) active pixels and (m, N + 1) where m = 1, . . . , M

1. A method of reconstructing a phase of a radiation wave fieldcomprising: determining a representative measure of an intensityvariation of the radiation wave field on a selected surface extendingglobally from one end to another of the radiation wave field;determining a representative measure of the intensity of the radiationwave field on the selected surface; transforming the representativemeasure of the intensity variation to produce a first representation ofintegral transform and to apply to the first representation of integraltransform a first filter and producing a first representation ofmodified integral transform; applying a first function of the firstrepresentation of integral transform to the first representation ofmodified integral transform and producing an untransformedrepresentation; applying a correction based on the measure of intensityon the selected surface to the untransformed representation;transforming the untransformed corrected representation to produce asecond representation of integral transform and applying a second filterto the second representation of integral transform; and applying asecond function of the second integral transform to the second modifiedrepresentation of integral transform to produce a measure of phase ofthe radiation wave field on the selected surface extending globally fromone end to another of the radiation wave field, wherein the first andsecond filters are measures of a representation of a Green's function ina space of eigenfunctions of the Helmholtz equation with separablecoordinates such that the first and second integral transforms areone-dimensional integral transforms in the space of the radiation wavefield.
 2. The method according to claim 1, wherein the first and secondintegral transforms are produced by using a Fourier transform.
 3. Themethod according to claim 2, wherein the Fourier transform is a fastFourier transform.
 4. The method according to claim 1, wherein the firstfunction uses a trigonometric function of the cosine type, and thesecond function uses a trigonometric function of the sine type.
 5. Themethod according to claim 1, wherein the first function uses atrigonometric function of the sine type, and the second function uses atrigonometric function of the cosine type.
 6. The method according toclaim 1, wherein the first and second filters are substantially thesame.
 7. The method according to claim 1, wherein the surfaces aresensibly planar.
 8. A computer program stored on a non-transitorycomputer readable medium that carries out the steps according toclaim
 1. 9. A computer program stored on a non-transitory computerreadable medium readable by a computer, including a means to carry outthe steps according to claim
 1. 10. An apparatus for reconstructing aphase of a radiation wave field comprising: a means for acquisition of arepresentative measure of an intensity variation of the radiation wavefield on a selected surface extending globally from one end to anotherof the radiation wave field; a means for acquisition of a representativemeasure of the intensity of the radiation wave field on the selectedsurface; a processor to sequentially perform the following steps: (i)transform the representative measure of the intensity variation toproduce a first representation of integral transform and apply to thefirst representation of integral transform a first filter and produce afirst representation of modified integral transform; (ii) apply a firstfunction of the first representation of integral transform to the firstrepresentation of modified integral transform to produce anuntransformed representation; (iii) apply a correction based on themeasure of intensity on the selected surface to the untransformedrepresentation; (iv) transform the untransformed correctedrepresentation to produce a second representation of integral transformand to apply to the second representation of integral transform a secondfilter; and (v) apply a second function of the second integral transformto the second modified representation of integral transform to produce ameasure of phase of the radiation wave field on the selected surfaceextending globally from one end to another of the radiation wave field,wherein the first and second filters are measure of a representation ofa Green's function in a space of eigenfunctions of the Helmholtzequation with separable coordinates such that the first and secondintegral transforms are one-dimensional integral transforms in the spaceof the radiation wave field.
 11. The apparatus according to claim 10,wherein the means of acquisition of a representative measure of theintensity variation of the radiation wave field includes at least amobile platform.
 12. The apparatus according to claim 10, wherein themeans of acquisition of a representative measure of the intensityvariation of the radiation wave field includes at least a lens ofvariable focal distance.